Lab Session 8

retaildata <- readxl::read_excel("retail.xlsx", skip = 1)
mytimeseries <- ts(retaildata[["A3349873A"]],
                  frequency=12, start=c(1982,4))
fit <- ets(mytimeseries)
summary(fit)
## ETS(M,A,M)
##
## Call:
##  ets(y = mytimeseries)
##
##   Smoothing parameters:
##     alpha = 0.5067
##     beta  = 1e-04
##     gamma = 0.3222
##
##   Initial states:
##     l = 63.0202
##     b = 0.7934
##     s=0.9391 0.912 0.9294 1.5281 1.0577 0.9868
##            0.9604 0.941 0.9431 0.901 0.9661 0.9353
##
##   sigma:  0.0488
##
##      AIC     AICc      BIC
## 4018.808 4020.494 4085.835
##
## Training set error measures:
##                     ME    RMSE     MAE         MPE     MAPE      MASE
## Training set 0.1576127 13.5932 8.70378 -0.09176148 3.711149 0.4596804
##                     ACF1
## Training set 0.008030422
autoplot(forecast(fit))

x1 <- window(mytimeseries, end=c(2010,12))
x2 <- window(mytimeseries, start=2011)
f1 <- snaive(x1, h=length(x2))
f2 <- hw(x1, h=length(x2), seasonal='multi')
f3 <- forecast(ets(x1), h=length(x2))
accuracy(f1,x2)
##                     ME      RMSE      MAE       MPE      MAPE     MASE
## Training set  7.772973  20.24576 15.95676  4.702754  8.109777 1.000000
## Test set     81.744444 100.00869 82.06667 20.549055 20.669738 5.143067
##                   ACF1 Theil's U
## Training set 0.7385090        NA
## Test set     0.6830879   1.67023
accuracy(f2,x2)
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set  0.2335309  8.817788  6.393003  0.1048658  3.20614 0.4006455
## Test set     78.8187005 95.779065 78.818701 20.0382885 20.03829 4.9395188
##                    ACF1 Theil's U
## Training set 0.06611033        NA
## Test set     0.52771581  1.630125
accuracy(f3,x2)
##                      ME       RMSE       MAE        MPE      MAPE
## Training set  0.3404863   8.784116  6.289258  0.1876561  3.128791
## Test set     95.8058083 113.202964 95.805808 24.5580948 24.558095
##                   MASE       ACF1 Theil's U
## Training set 0.3941439 0.01203982        NA
## Test set     6.0040903 0.61158849  1.922081
bicoal %>% ets %>% forecast %>% autoplot

chicken %>% ets %>% forecast %>% autoplot

dole %>% ets %>% forecast %>% autoplot

usdeaths %>% ets %>% forecast %>% autoplot

bricksq %>% ets %>% forecast %>% autoplot

lynx %>% ets %>% forecast %>% autoplot

ibmclose %>% ets %>% forecast %>% autoplot

eggs %>% ets %>% forecast %>% autoplot

ausbeer %>% ets %>% forecast %>% autoplot

Lab Session 9

autoplot(usnetelec)

No transformation required

autoplot(mcopper)

(lambda <- BoxCox.lambda(mcopper))
## [1] 0.1919047
mcopper %>% BoxCox(lambda=lambda) %>% autoplot

autoplot(enplanements)

(lambda <- BoxCox.lambda(enplanements))
## [1] -0.2269461
# I don't like such strong transformations. Will use 0 instead
enplanements %>% BoxCox(lambda=0) %>% autoplot

autoplot(a10)

(lambda <- BoxCox.lambda(a10))
## [1] 0.1313326
a10 %>% BoxCox(lambda=lambda) %>% autoplot

autoplot(cangas)

(lambda <- BoxCox.lambda(mytimeseries))
## [1] 0.1276369
mytimeseries %>% BoxCox(lambda=lambda) %>% autoplot

f4 <- stlf(x1, lambda=lambda, h=length(x2))
accuracy(f4,x2)
##                     ME      RMSE       MAE        MPE      MAPE     MASE
## Training set -0.604219  7.214593  5.176787 -0.2755991  2.525501 0.324426
## Test set     79.914332 97.176605 79.974531 20.2410238 20.263570 5.011954
##                    ACF1 Theil's U
## Training set 0.02490237        NA
## Test set     0.57946539  1.643056
autoplot(f4) + autolayer(x2, series="Test data")

Lab Session 10

usnetelec %>% autoplot

usnetelec %>% diff %>% autoplot

usgdp %>% autoplot

usgdp %>% diff %>% autoplot

lambda <- BoxCox.lambda(mcopper)
mcopper %>% BoxCox(lambda=lambda) %>% autoplot

mcopper %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot

enplanements %>% log %>% autoplot

enplanements %>% log %>% diff(lag=12) %>% autoplot

enplanements %>% log %>% diff(lag=12) %>% diff %>% autoplot

visitors %>% autoplot

lambda <- BoxCox.lambda(visitors)
visitors %>% BoxCox(lambda=lambda) %>% autoplot

visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% autoplot

visitors %>% BoxCox(lambda=lambda) %>% diff(lag=12) %>% diff %>% autoplot

mytimeseries %>%
  BoxCox(lambda=0.12) %>%
  diff(lag=12) %>%
  diff(lag=1) %>%
  autoplot

Lab Session 11

wmurders %>% autoplot

wmurders %>% log %>% autoplot

fit <- auto.arima(wmurders, lambda=0, stepwise=FALSE, approximation=FALSE, max.p=10, max.q=10, max.order=10)
checkresiduals(fit)

##
##  Ljung-Box test
##
## data:  Residuals from ARIMA(1,2,1)
## Q* = 10.562, df = 8, p-value = 0.2278
##
## Model df: 2.   Total lags used: 10
forecast(fit) %>% autoplot

wmurders %>% ets %>% forecast %>% autoplot

wmurders %>% ets(lambda=0, model="AAN") %>% forecast %>% autoplot

Lab Session 12

lambda <- BoxCox.lambda(mytimeseries)
auto.arima(mytimeseries, lambda=lambda)
## Series: mytimeseries
## ARIMA(1,1,1)(0,0,2)[12] with drift
## Box Cox transformation: lambda= 0.1276369
##
## Coefficients:
##          ar1      ma1    sma1    sma2   drift
##       0.5676  -0.9565  0.9543  0.5450  0.0096
## s.e.  0.0523   0.0183  0.0528  0.0357  0.0022
##
## sigma^2 estimated as 0.02727:  log likelihood=140.48
## AIC=-268.95   AICc=-268.73   BIC=-245.31
(arimamod <- auto.arima(mytimeseries,
              D=1, lambda=lambda,
              stepwise=FALSE,
              approximation=FALSE))
## Series: mytimeseries
## ARIMA(1,1,2)(0,1,1)[12]
## Box Cox transformation: lambda= 0.1276369
##
## Coefficients:
##           ar1     ma1      ma2     sma1
##       -0.4147  0.0703  -0.3325  -0.5811
## s.e.   0.1898  0.1822   0.0725   0.0555
##
## sigma^2 estimated as 0.009902:  log likelihood=326.42
## AIC=-642.84   AICc=-642.67   BIC=-623.3
newretaildata <- readxl::read_excel("8501011.xls", sheet=2, skip=9)
mynewdata <- ts(newretaildata[["A3349873A"]],
                  frequency=12, start=c(1982,4))
mynewdata <- window(mynewdata, start=2014)

lambda <- BoxCox.lambda(mytimeseries)
(etsmod <- ets(mytimeseries))
## ETS(M,A,M)
##
## Call:
##  ets(y = mytimeseries)
##
##   Smoothing parameters:
##     alpha = 0.5067
##     beta  = 1e-04
##     gamma = 0.3222
##
##   Initial states:
##     l = 63.0202
##     b = 0.7934
##     s=0.9391 0.912 0.9294 1.5281 1.0577 0.9868
##            0.9604 0.941 0.9431 0.901 0.9661 0.9353
##
##   sigma:  0.0488
##
##      AIC     AICc      BIC
## 4018.808 4020.494 4085.835
f1 <- snaive(mytimeseries, h=length(mynewdata))
f2 <- hw(mytimeseries, h=length(mynewdata), seasonal='multi')
f3 <- forecast(etsmod, h=length(mynewdata))
f4 <- stlf(mytimeseries, lambda=lambda, h=length(mynewdata))
f5 <- forecast(arimamod, h=length(mynewdata))
checkresiduals(f1)

##
##  Ljung-Box test
##
## data:  Residuals from Seasonal naive method
## Q* = 804.66, df = 24, p-value < 2.2e-16
##
## Model df: 0.   Total lags used: 24
checkresiduals(f2)

##
##  Ljung-Box test
##
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 40.245, df = 8, p-value = 2.884e-06
##
## Model df: 16.   Total lags used: 24
checkresiduals(f3)

##
##  Ljung-Box test
##
## data:  Residuals from ETS(M,A,M)
## Q* = 26.597, df = 8, p-value = 0.0008296
##
## Model df: 16.   Total lags used: 24
checkresiduals(f4)

##
##  Ljung-Box test
##
## data:  Residuals from STL +  ETS(M,A,N)
## Q* = 69.123, df = 20, p-value = 2.531e-07
##
## Model df: 4.   Total lags used: 24
checkresiduals(f5)

##
##  Ljung-Box test
##
## data:  Residuals from ARIMA(1,1,2)(0,1,1)[12]
## Q* = 16.502, df = 20, p-value = 0.685
##
## Model df: 4.   Total lags used: 24
accuracy(f1,mynewdata)["Test set","RMSE"]
## [1] 77.30981
accuracy(f2,mynewdata)["Test set","RMSE"]
## [1] 52.98594
accuracy(f3,mynewdata)["Test set","RMSE"]
## [1] 53.71773
accuracy(f4,mynewdata)["Test set","RMSE"]
## [1] 39.35806
accuracy(f5,mynewdata)["Test set","RMSE"]
## [1] 29.27015
autoplot(f5) +
  autolayer(mynewdata, series="New data")